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qq , We study the statistical mechanics of classical and quantum systems in non-equilibrium steady states. Emphasis 

is placed on systems in strong thermal gradients. Various measures and functional forms of observables are presented. 
The quantum problem is set up using random matrix techniques, which allows for the construction of the master 
■ equation. Special solutions arc discussed. 
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I. INTRODUCTION 
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The efforts to develop statistical mechanics of non-equilibrium systems are easily traced back to the foundations of 
the subject. As statistical mechanics was developed, it was a natural step to consider how to extend the theory. One 
can see this in the handwritten notes of J.W. Gibbs |5J], in which he formulated the theory of statistical mechanics. 
Gibbs certainly expended some effort to understand how he might characterize the non-equilibrium steady-state. 
Gibbs certainly puzzled over how to define entropy for systems in non-equilibrium steady states, as well as how to 
define the statistical mechanics of a gas in a pressure or temperature gradient. However, no solutions were found. 
The field has remained active since that time. In the past ten to fifteen years, the emergence of results in classical 
chaos has led to renewed interest and many new ideas have come about || . 

Non-equilibrium statistical mechanics is replete with unanswered questions. While many theoretical techniques have 
been suggested to treat transport, a general understanding of the statistical mechanics is still lacking. In these lectures 
we survey some recent results on non-equilibrium steady states in classical and quantum systems using dynamical 
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II. SOME RESULTS FROM CLASSICAL NON-EQUILIBRIUM STATISTICAL MECHANICS 

We examine model systems in which a Hamiltonian is coupled to heat reservoirs at two endpoints, as shown in 
Fig. |. 
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FIG. 1. Model system of a Hamiltonian H coupled to two thermal reservoirs, at temperatures T co id 
separated by a distance L in the a;— direction. 



T? and T hot = T 2 ° 



The dynamics will be that of the given Hamiltonian, except at the boundaries where it is coupled to the reservoirs. 
The boundary conditions can be implemented in a variety of ways S. Before we discuss quantum systems, we present 
some results for classical systems, which summarizes results in Refs. fj- Q and references there in. 

Thermal Conductivity. 

The transport coefficient in this case is the thermal conductivity n. There are two approaches that are used to 
compute k. The first is the Green-Kubo approach, where the non-equilibrium transport is computed in terms of 
equilibrium auto-correlation functions. In this case: 

K (T) = ±Jllt Jdx (J(x,t)J(x ,0)) BQ . (1) 

The heat flow J can be readily computed from the energy-momentum tensor, J — T 01 , for heat flow in the x— direction. 
Alternately, one can compute the conductivity directly in a non-equilibrium steady state using Fourier's Law. Mea- 
suring the local heat flow J inside the system, and the local temperature T(x), leads to 

<T) = _(Jh™ t (2) 

V ' WT(x) V ; 

For many systems k(T) behaves as a power law, or can be approximated by a power law in temperature ranges of 
interest: 

= ^. (3) 

We should make a few notes concerning Fourier's law. 

• A constant gradient inside a system does not guarantee that Fourier's law holds. When 7 is very small, as in 
the Td Fermi-Pasta- Ulam (FPU) 0— model around T ~ 1, the system can already be too far from equilibrium 
for Fourier's law to hold, yet the nearly constant k still provides a near constant gradient and the illusion that 
the law holds. In the case of the FPU model one finds that Fourier's law is not valid even locally when the 
system has a near constant temperature gradient. 

• A strongly-curved temperature profile T(x) (as in Fig. 2) does not mean that Fourier's law is not valid. It has 
been shown that even strongly curved profiles can be derived analytically if Fourier's law is satisfied locally. 
So to contrast with the previous point, the curvature of T(x) does in itself provide no information on whether 
Fourier's law is valid. 

• A system does not need to have a bulk limit to satisfy Fourier's law. In the 1-d FPU model, there is no bulk 
limit, and n depends on the system size L. Never the less, Fourier's law is valid even for reasonably strong 
thermal gradients. 

• Fourier's law (2) and the Green-Kubo result (I) have been seen to agree and describe the properties of systems 
even very far from equilibrium, for cases where T(x) is strongly curved. Hence linear response theory predictions 
are quite robust and can be expected to hold in strong non-equilibrium environments. 
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Temperature Profile. 



When we discuss the temperature profile inside, T(x), we must distinguish the temperatures of the heat reservoirs, 
denoted and T§ , and the temperatures just inside the system T\ and This is indicated in Fig. ^. The boundary 
jumps ST can be quantitatively understood, as we mention below. The temperature inside T(x) is given most directly 
in terms of the extrapolated temperatures T\^- 




FIG. 2. Temperature profile of a system such as that depicted in Fig. 1. The boundary conditions impose temperatures 
T co i d = Tf and T hot = T 2 ° on the ed ges, x = 0, L. Just inside the system, there is a boundary jump ST. Extrapolating the 
smooth temperature profile T(x) to the edges defines the temperatures Ti and Ti. 



We can understand the curvature in the temperature profile using Fourier's law. When the thermal conductivity 
can be described by a power law in the temperature range of interest (Ti, T2), such as k(T) = AT~ 7 , we can integrate 
Fourier's law and re-express the result in terms of extrapolated temperatures (T 1; T 2 ): 



T(x) = 




1-7 



7^1 
7-1 



(4) 



Note that k(T) need not have a bulk behavior, nor does it need to behave globally as a power law to describe T(x) 
in non-equilibrium systems. These formulas are quite robust. 



Boundary Jumps. 

An interesting and often overlooked property of non-equilibrium steady states is the presence of boundary jumps. 
In fluids which are sheared by moving walls, one observes that the velocity just inside the fluid can be different from 
the velocity of the wall. For systems in a thermal gradient, there is an analogous jump in the value of the temperature 
just inside the system. This jump, ST, is illustrated in Fig. 2. 

The jump arises due to a mismatch of the mean free path with the edge of the system. It has the general form 



n dT 
STi = — Ti — rj — — 
on 



(5) 



boundary 



where 77 is a parameter on the order of the mean free path I and n denotes the normal to the boundary |Tc| ]. This 
relation has been verified in two independent manners. 

In the first, the boundary jumps ST have been shown to be linearly related to the normal derivative of the tem- 
perature profile T(x), extrapolated to the edges of the system. The coefficient 77 can then be extracted. This relation 
has been verified in several systems. In the (/> 4 model in 1-d, for example, 



n(T) = (6.1±0.5)T 



-1.5±0.1 



(0) 



To check if we have the correct understanding, we can make an estimate of 77 from kinetic theory. Recall that the 
thermal conductivity is related to the mean free path by k ~ Cyct by a standard kinetic theory argument. For 
systems where Cy and c are largely temperature independent, such as in various lattice models, we expect k ~ I. 
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Dimensionally, the coefficient r\ should be on order of the mean free path £, so that r\ ~ k. Both the magnitude and 
power law behavior are consistent with the thermal conductivity obtained in the (j) A model in 1-3 D as well as the 
FPU model. 

An independent verification of the behavior of the jumps can be made by studying how the jumps depend on the 
heat flux. We let 77 = qk, where a is a constant to be determined. From Eq. (^]), we can then associate the heat flux 
with the right side of Eq. (j^) : 

T< - TP ~ a(J) NE ~ a(T 2 ° - + ■■■. (7) 

a is a coefficient which measures the efficacy of the boundary conditions. The linear relation between ST and (J)ne 
has been verified, a has been extracted for a variety of boundary thermostats, and can show strong dependence on 
those boundary conditions. The understanding of these jumps together with the temperature profile ([|) provides a 
complete description of T(x) in terms of the boundary temperatures, (T-P, T® ). 



Observables. 



If we probe an observable somewhere in the system in the non-equilibriums steady state, can we understand 
generically how we expect these to behave? It is a natural idea that a physical observable A will deviate from its 
value in local equilibrium as we move further away from equilibrium. In our case, a temperature gradient VT provides 
a measure of how far we are from equilibrium. The deviation of an observable SA from its local equilibrium value is 
expected to behave as 

This expansion is in even powers since the result should be insensitive to which side of the box is taken as the hot 
side and which is cold. (Using Fourier's law, we can equally well expand in terms of the current J.) While this seems 
natural, there has been some contention on analogous expansions in sheared fluids. There non-analytic dependencies 
on the shear rate have been observed in some numerical simulations, but the situation has not been entirely clarified 
jnj. The expansion coefficients C.4, C^, . . . are in principle dependent on physical quantities, such as T and L, the 
size of the lattice in the direction of the gradient. If the departures from local equilibrium are local, then the the 
coefficients would be expected to be independent of L. It has been shown that even in systems with bulk transport 
behavior, the coefficients do indeed depend on the system size. 



Far from Equilibrium Spatial Dependence. 

The expansion in (8) has the advantage that we can combine it with (4) and obtain the spatial information on how 
observables vary from local equilibrium values within a system: 

SA „ ( J \ 2 „ ( 1 

where 



A =C +\jf J =^U+fa» (9) 



j>1-7 A 

a =777 , 6 = 7-1. (10) 



(J) 



NE 



So knowing C4, we can also predict the spatial variation of the non-equilibrium observable A. This has been tested 
and verified in classical lattice models in d = 1 — 3 dimensions. 

The form of (9) brings to light an interesting relation to coarse graining. Let us assume we have an observable 
which is locally non-equilibrium, so that SA/ A at some point x' is non-zero. Using traditional arguments, we would 
take a larger box containing x' and assume that on the larger scale, the behavior is closer to local equilibrium. In 
examining (9) we see that it is a positive definite function of position (up to the overall sign of C4). So any averaging 
over a larger box will not necessarily converge to SA/ A = 0, required for local equilibrium. 



Breaking of Local Equilibrium. 



Local equilibrium is a property which is usually invoked to allow application of equilibrium physics to a problem of 
interest. It is seldom quantified. Using the expansion (8), it is possible to quantify the 'quality' of the local equilibrium 
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assumption. The local temperature at position x in the model of Fig. 1 is given by Tk = (py- A natural measure for 
the deviations from local equilibrium is the deviation of the momentum distribution from the Maxwellian distribution. 
Consider 

A=(p 4 ). (11) 

In local equilibrium, we expect a Gaussian distribution of momenta, thermalized at the local temperature T(x) (or 
Tk), so that 

A eq = (p 4 ) eq = 3(p 2 ) 2 = 3T(x) 2 . (12) 

Hence a measure of local equilibrium is given by 

S(pi} = (pi)-3(pi) 2 = ((pi)), (13) 

where ((pf )) is the cumulant at location k. In principle, we can examine the behavior of other cumulants are measures 
of the breaking of local equilibrium, such as ((pf)) = (pf) — lb{p\){p\) + 30(pf) 3 . It is convenient to normalize the 
cumulants by the local temperature, ((p"))/T™ //2 . These provide a quantitative measure on how far we are from local 
equilibrium. In local equilibrium, ((p n )) =0 (n > 2). 

Here a few properties which have been observed in lattice models: 

• Contrary to naive intuition, the steepest gradient in Fig. 2 does not lead to the system being furthest from 
local equilibrium. In fact, the converse is true — the system is furthest from equilibrium in the flattest region. 

• It has been found that VT/T provides a good measure of how far we are from equilibrium and that the cumulants 
behave as 

In the <f> 4 and FPU models,^ = l.l(8)i°- 9 < 2 ) (T = 1) and C™ = 4.3(4)L°-"( 2 ) (T = 8.8). There seems to 
be a weak T dependence for Cle which is difficult to establish. 

• These results are consistent with d > 1 in lattice models at the same temperatures. 

• Using C^ E , C^ u , one can predict the shape of ((p 4 }}/T 2 in the system using (9), which has been verified 
numerically. 

• As the system moves away from equilibrium by increasing the difference in the boundary temperatures, each 
point in the interior deviates from local equilibrium in a predictable manner, without any threshold. Away from 
equilibrium, local equilibrium is an approximation that is quite good for small gradients since the deviations 
from it only vary as (VT) 2 . 

• Similar results hold for higher momentum cumulants. 

There are many open questions here that related to measurements in lattice models. The naive expectation would 
be that since the gradients and the cumulants are local, their relationship should not depend on the system size L. 
There does seem to be size dependence, even in systems with bulk transport, many mean free paths away from the 
boundaries. 



Corrections to Linear Response. 

There is an interesting relation between linear response predictions for transport coefficients and the notion of 
local equilibrium. For instance, if there is a departure of a predicted transport coefficient from a measured value for 
a system far from equilibrium, is this an indication that (a) corrections to linear response are needed or (b) local 
equilibrium is no longer a valid assumption. If the latter is true, temperature can no longer be unambiguously defined. 
Near global equilibrium, Fourier's law is satisfied globally so that there is a constant thermal gradient: 

Jo - -k(T)(T 2 - Ti)/L ~ -k(T)(T 2 ° - T?)/L. (15) 
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As the temperature difference increases, curvature in T(x) usually develops, and Fourier's law is satisfied locally. For 
k <~ T~ 7 , Fourier's law can be integrated to obtain the next leading order correction to Jo, due to curvature in T(x). 
Denote this heat flow, Jlr, which includes the curvature correction. It behaves as: 



Jlr = Jo 



1 + 



7(7- 



24 



VT 



+ 



(16) 



Notice that the curvature correction to Jo behaves as L 2 . 

Using numerical simulations on lattice models, one can push the system very far from equilibrium and examine how 
J behaves. It is found that as the gradient increases, the energy that can be pumped through the system becomes 
less than that predicted by linear response theory even when it is applied locally. The 'violation' of linear response is 
defined as Slr and is found to behave as 



Slr = 



Jlr 



\ T I 



(17) 



In the (ft 4 and FPU models, it is found that C 



lr 



4(3) L i.o(2) ( T = ^ 



and C[ PU 



-6.6(8)L - 9 « (T = 8.8). 



Notice that this will naively give rise to decreasing current J with increasing gradient. However, when the gradient 
is this large, the higher order terms in (VT/T) becomes as important. It is believed that the current will eventually 
saturate under extremely large gradients. 

A very interesting picture emerges. Within error, the violation of linear response and local equilibrium are closely 
connected, occurring in the same manner: 



OLE 



JLR- 



(18) 



Further, local equilibrium and linear response have no threshold, and break down in a similar manner as the system 
moves away from global thermal equilibrium. 



Non-Equilibrium Equation of State. 

For a simple one-component theory, the equation of state will have a simple form such as P eq {T) or E eq (T), where 
P and E are the equilibrium pressure and energy density. In the non-equilibrium steady state, one expects new 
variables to emerge. For thermal gradients, VT will become an independent variable. What is surprising, is that 
even in systems with bulk behavior, there can develop a system size dependence for the non-equilibrium 'equation of 
state'. For (ft 4 and FPU lattice models, it is found that 

P(T,VT,L) = P eq (T) 

where Cf, = 1.5(1.2)L°- 9 W (T = 1), C FPU = 4.1(6)T°- 30 ( 4 ) (T = 8.8). Notice that the L dependence makes the 
equation of state non-local even though it is measured locally, in a system with bulk behavior. If we study the energy 
density, we have a similar result: 

E(T, VT, L) = E eq (T) 

where C% = 0.5(3)T°' 9 ( 2 ) (T = 1), C£ PC/ = 1.7(7)T°' 3 W (T = 8.8). In contrast, Extended Irreversible Thermody- 
namics (EIT), predicts a local behavior for P — P eq in contrast to our observations. 



1 + Cp 



VT 



(19) 



1 + C E 



/VT 



(20) 



Dimensional Loss. 

A consistent picture that emerges from classical simulations of non-equilibrium steady states is that the accessible 
volume of phase space contracts onto a fractal set. The dimensionality is obtained from the Kaplan- Yorkc estimate, 
which involves calculating the full spectrum of Lyapunov exponents. There is a corresponding loss of dimension AD. 
An important question that arises in non-equilibrium systems is whether the dimensional loss is extensive, and persists 
in the bulk, or simply a manifestation of low dimensionality. In systems with bulk behavior, we expect AD to behave 
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"extensively" . That is, AD should remain relatively the same under the same local non-equilibrium conditions when 
we change the size of the system. Extensivity has been shown in the following form: 

^-=C D J* (21) 

where Vi n is the volume of the system inside, excluding the thermostatted degrees of freedom, while D includes all 
degrees of freedom. The coefficient is given by 

Cd = — r^g — =2(1 + -^— - ) m eq 2 (22) 

Here \%\ ax is the maximum Lyapunov exponent of the full system plus bath, in thermal equilibrium. Some results are 

• Extensivity is explicitly related to the macroscopic transport properties through entropy production, which leads 
to (22). The first systematic study of the role of dynamical thermostats on dimensional loss was also shown in 
[8]. 

• While it was emphasized that the extensivity of AD is not compatible with local equilibrium so that it is 
questionable jL2|. As was mentioned above, violations of local equilibrium appear in the same manner as 
dimensional loss. Hence there is no conflict. 

• An estimate of dimensional loss in an ideal gas gives an idea of the magnitude of the effect. Using the standard 
estimates of \ ma x 0] yields Cd — [9pv 2 ln(4Z/d)/2] _1 , where p is the density, v is the average particle velocity, I 
is the mean free path and d is the particle diameter. For VT/T ~ 0.01 m _1 , AD ~ 10 9 m~ 3 at room temperature 
— quite large, yet far smaller than the total number of degrees of freedom. 

• The results are satisfying from the physics point of view. AD pertains to the whole system, it includes the 
temperature profile which is curved in general, boundary temperature jumps and the various types of ther- 
mostats. Yet, AD can be related to macroscopic transport with the thermostat dependence cleanly separated 
into X m ax- Furthermore, AD has extensive behavior with respect to the internal volume wherein the the system 
is manifestly in non-equilibrium. We have seen that AD and \%\ ax both depend strongly on the thermostat 
used, while their product ADA^ aa . is thermostat independent and related to macroscopic quantities. 

• ^max i s n °t unique: In global thermal equilibrium, different choices of heat-baths can lead to very different 
values. The result is that dimensional loss is not unique either, only the product behaves macroscopically and 
can be related to thermodynamic quantities 



Non-Equilibrium Temperature Rcnormalization. 

In expansions such as (8), it is natural to consider whether other powers are important, such as (V n T)/T. In the 
region where the temperature dependence of the thermal conductivity can be described by a power law, k = AT 1 , 
one can show, using Fourier's law, that 

^<^ ( „VnV^Y (23) 

fe=i v J 

While Fourier's law holds only close to local equilibrium, the deviations from it is of order (VT/T) 2 so that the 
difference is a higher order correction in the expansion in (VT/T). Hence the expansion in (8) is sufficiently general. 

When local equilibrium is broken, the definition of T is no longer unique. For the observables discussed above, it is 
possible see that certain redefinitions of T leave the observables invariant to leading order. For the generic redefinition 



T = T' + v\—j , C' A = C A + ^-), (24) 

the non-equilibrium deviations of J, P, E are affected covariantly. Further, the local equilibrium violations seen in 
((p"))/T™/ 2 are invariant under such redefinitions, up to the order we consider. The physics, of course, is invariant 
under any redefinition in temperature. 

Summary of Classical Results. 
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We have seen that we can begin to understand the behavior of systems both near and far from equilibrium. 
Certainly many questions remain, including the behavior of phase transitions in strongly non-equilibrium steady-state 
environments, which seem to display interesting differences with the equilibrium understanding of phase transitions [7] 
. But we would now like to turn to the treatment of the quantum analogs of the above results. 

III. DYNAMICAL QUANTUM HEAT BATHS 

Before we explore the quantum analog of Fig. 1 and the results for classical non-equilibrium steady states, let us 
first consider how to formulate a Hamiltonian description of a heat bath. We would like to use chaos in a similar 
manner; the (chaotic) interaction with the baths should cause the system of interest to thermalize. We decompose 
the system+heat bath problem as: 

H(XP,xp) = H S (XP) + H^X-xp) (25) 

= H s [XP) + H B (xp) + W(X- 1 xpj . 

The system, Hs(XP), and the bath, Hb(xp), Hamiltonians are coupled through the interaction W(X;xp). In this 
formulation, for simplicity, the interaction does not include any system momentum dependence. The density of states 
of the bath+interaction is defined by Hi through 

p(E) = TWffli - E). (26) 

In order to use chaotic methods, we should make a few remarks concerning Random Matrix Theory (RMT), which 
has arisen as the quantum counterpart to classical chaos. 

What is a random matrix? A matrix with random matrix elements JD^O ], often taken as Gaussian random numbers. 
This is convenience rather than necessity, as the Gaussian model is far more tractable and well developed than other 
choices. The motivation is simple. The Schrodinger equation, H\ipi) — E i \^i i ) 1 cannot be solved for complex many- 
body systems. This lead Wigner (1951) to adopt a stochastic description instead of a dynamical description and to 
replace the Hamilton operator by a random operator. 

How are the matrix elements of a random matrix Hij distributed? Assume that 

• H is invariant under "rotation" (no preferred basis in Hilbert space), 

• the matrix elements H+j (l<j<N) are independent. 

It can then be shown that H^ are Gaussian random variables. Depending on the symmetries of the problem, there are 
only three possible "rotations" in Hilbert space: orthogonal, unitary and symplectic transformations. Accordingly, 
there are only three Gaussian random matrix ensembles denoted GOE, GUE and GSE (note that generalizations exist, 
for example chiral ensembles in QCD). For GOE, P({Hij}) — K exp[— - Hfj/ ' (HijHji)]. In general the variance 
(HijHji) depends on 

The classical analogy is the velocity distribution in the kinetic theory of gases. Newton's equations describing the 
(large number of) molecules of an ideal gas cannot be solved. This motivated Maxwell (1859) to use, for the first time 
in physics, a statistical approach to treat a dynamical problem (— ► random vector). He assumed that 

• space is isotropic (no preferred direction in velocity space), 

• the components Vi (i — 1, 2, 3) of the velocity vector are independent. 

He then found that the velocities are Gaussian distributed (the well-known Maxwell distribution), P({vi}) = 
Kexp{- vf/(v?)}, with (vf) = 2kT/M. 

Heat Bath Model. 

We will use the eigenbasis of Hb(xp) to define the matrix elements of W . If we take the interaction to be chaotic, 
it is reasonable to choose the statistical ensemble to be defined through W according to the GOE ensemble. As the 
Gaussian ensemble is entirely characterized by its first two moments, this is taken to be in the form: 

[W(X)] kl = 0, [W(X)] k i[W(Y)] mn = [S km 5 nl + S kn S lrn ]g kl (X - Y). (27) 
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The physics is put into the properties of the correlations. We use the form: 



Sa(X) = 



V£G(X/X ). 



(28) 
(29) 



This type of parameterization was introduced some time ago in the study of heavy ion collisions p5| and widely used 
afterwards [|l6-19|. The physical motivation is: 

• The interaction is expected to couple only neighboring energy states — ► W = random band-matrix with band- 
width kq <C N. For an ordinary random matrix (not banded), Q a b 2 = constant. 

• With increasing energy, the states \a) and \b) become more and more complex. Consequently, their overlap is 
expected to decrease. This is the physical origin of the factor containing the density of states of the bath. 

• It is possible that correlations exists between W at different values of position X . This information is included 
in the correlation function G(X/Xq), where Xq is the characteristic length scale. 

We emphasize that the coupling to the random matrix bath is fully determined by the variance. The parameters that 
enter are: - the spreading width; Xq - correlation length measuring how rapidly the system decorrelates due to 
the bath; kq - band width of the random matrix W; p = po exp[/3e] - density of states of the bath used to define the 
thermodynamic temperature through (3 = l/T = dlogp/de. 



Influence Functional Approach. 

The interaction of the system Hs with a chaotic bath Hb + W can act to thermalize the system. The evolution 
equation for the system can be obtained by integrating over the bath. This is most readily achieved using the influence 
functional [^0| . The density matrix of the system+bath evolves according to 

p(t)=J{t,t')p{t'). 

The evolution operator J has a simple path integral representation which is the squared amplitude: 
J(X, x, X', x', t\Y, y, Y', y', t 1 ) = K(X, x, t\Y, y, 0)K*(X', x', t'\Y', y', 0) 

Dx 



' Y 

Dx' 



DXexp 

x' 



-S(X,x) 



(30) 

(31) 
(32) 



y' 



DX' exp 



■-S(X',x') 



S is the action for the system+bath. The density matrix can be expressed in terms of this path integral averaged 
over the bath degrees of freedom p0[ 



p(X,X',t)= / dX dY po(Xo,X' ) 



X(t)=X 



X(0)=X o 



VX(t) 



X'(t)=X' 



X'(0)=X^ 



VX'(t) 



x exp 



[So(X(t))-S (X'(t))} 



£(X(t),X'(t),t). 



with the influence functional given by: 

C(X(t),X'(t),t) - li\ (T a exp 



dt" H^X 1 {t")) 



Texp 



dt'HxiXit')) 



(33) 



(34) 



Here T (T a ) is the (anti) time ordering operator. At this point we can expand the propagators in the influence 
function, and use the statistical properties of the matrix elements. This allows a systematic approach to include finite 
temperature effects into the quantum dynamics. To first order in /? = l/T one finds 



£{X,X') = exp 









H 







x exp 



i(3T l 



' X( S )-X'( S ) 



ds[X {s)+X'{s)]G' 



1 



X(s)-X'(s) 
X 



(35) 
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where G'(x) — dG(x)/dx. From this one can derive the effective time evolution of the density matrix for the system: 

(Px-Px>) (36) 



+ iT l 



p(X,X',t), 



with an arbitrary initial condition p{X, X',0) = po(X,X'). This equation provides the quantum dynamics of the 
system at finite temperature p0| ]. The initial conditions can be taken as gaussian, ipa(X) = exp[— X 2 /4cr 2 ]/[27rcr 2 ] 1 / 4 , 
so that the initial density matrix is 



Po(X,X') 



1 



V27R7 2 



-(X 2 +X' 2 )/4a 2 



1 



-(4r 2 +s 2 )/8<r 2 



V27TCT 2 



where 



X +X' 
r = , s = X - X'. 



(37) 



(38) 



In order to learn about the transport behavior, it is convenient to define the characteristic function d(s, k, t), whose 
moments give the transport coefficients: 



p(r, s,t)= I — exp f — ) d(s, k, t). 



Then the cumulants are extracted as follows: 

1 fik 



n=l 



lnd( S ,M)| s=0 = E ^ ( j) ^d(s,k,t)\ k=0 = E^[ (£) n d pn 



(39) 



(40) 



n=l 



The first cumulant is just the average, ((r)) = (r), the second is ((r 2 )) = (r 2 )— (f*) 2 , the third is usually referred to as 
kurtosis, and so forth. The cumulants are closely related to transport coefficients. 



IV. FREE PARTICLE STRONGLY COUPLED TO THE RESERVOIR 



The master equation can be solved in many limits. To see the effects of the bath, we take the potential U(X) = 0. 
Then density matrix is expressed as: 



p(r,s) = e lkr / h y{s). 



This leads to the simplified master equation: 



ihk d . ( s \ d 



*(s) = 



which can be solved through direct integration. For a general correlation function G: 

. . f ikr MT l f s , , 1 - G 

p(r lS )=ex P |— --^J Q ds - x {lXoM/k)G , 

Consider some examples. 

Exponential Correlator: G(x) — exp(— \x\) Defining a = jXgM/k and b = MT^/hk 



*(*) 



1 + a 



-(l+a)b/a 



D sb/a 



s a + exp(s) 

Cosine Correlator: G(x) = cos(x) In this case, the solution is of the form 



(41) 



(42) 



(43) 



(44) 
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= ( l + (a+^^taW2)\ 
W ^1 + (a- V^^T tans/2)/ V 7 

From these explicit solutions to the strongly coupled master equation, one can compute the coordinate and momentum 
cumulants and demonstrate how the quantum system evolves in time in the presence of this type of bath. As our 
interest is in the weak coupling limit where the bath looks closer to an idealized heat bath, we will not pursue this 
here. 



V. RANDOM MATRIX THEORY MASTER EQUATION: 
WEAK COUPLING PLUS ONE BATH 



A. Hamiltonian 



We consider a system S coupled to a heat bath B. The corresponding Hamiltonian is given by 

H = H S + H B + W , (46) 

where Hs — P 2 /2M + U(X) describes the system [a particle moving in a potential U(x)], Hb describes the bath and 
W = X ® V is the coupling between the system and the bath. The coupling Hamiltonian W is taken linear in the 
position x of the system and V is an operator acting on the bath. We denote by e a the energy eigenstates of the bath, 
Hb\o) — s a \a), and assume that the heat bath is complex (chaotic). It is then justifiable to model V by a random 
matrix: For a complex system, the states \a) are expected to be highly complex and to display random features; it is 
thus reasonable to assume that V a b is a stochastic function of a and b. 



B. Master equation 



Our aim in this section is to provide an alternate derivation of a quantum master equation which gives the time 
evolution of the (ensemble averaged) reduced density operator of the system, ps(t) — tTBp(t), when the coupling is 
weak. Here p(t) is the total density operator for system plus bath. We shall employ a method well-known in quantum 
optics (see for instance Starting point of the derivation is the von Neumann equation for p(t) written in the 

Interaction Picture, 

^ = -i[W(t),p(t)}, (47) 
with A(t) = exp(iH t)Aexp(—iH a t) and H Q = Hs + Hb- Equation (pTjj) can be formally integrated to give 



p(t)=p(0)-i dt' [W(t'),p(t')] . (48) 
Jo 

We now substitute this expression for p(t) inside the commutator of ( ff7| ) and obtain 

^ = -i[W(t),m)- l^dt'[W(t),[W(t%m}} • (49) 

This equation is still exact. We next assume that initially the system and the bath are not correlated and that the 
latter is in thermal equilibrium, i.e., p(0) = p s (0) <g> p B (0) with p B {0) = Z^expi-PHg)- We also assume (Born 
approximation) that the bath remains in equilibrium at all times so that we can write p(t') — ps(t') ® Ab(0). Clearly, 
this approximation is true provided the coupling between the system and the bath is weak. We now take the trace 
over the bath and note that trspit) = exp(iHst)tTBp(t) exp(-iHst) = ps(t). After ensemble averaging we obtain 

^2 = _ J* drK{r){x{t)x{t - r)p s (t - r) x(t - r)p s (t r)x(t)} 

+ [ dTK(-r){p s (t - r)x(t - r)x{t) - x(t)p s {t - r)x(t - r)} , (50) 
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where we have used W(t) — x(t) ® V(t) and introduced the bath correlation function K(t) = K'[t) + iK"(t) 
trB[V(t)V(0)pB(0)]- Its ensemble average is given by the Fourier transform of V a b 2 with respect to e\, |2^| , 

" + OO 



(51) 



Transforming Eq. (^9|) back into the Schrodinger Picture yields 

dp s (t) _ 



dt 



i[H s ,p s (t)]- / dr[x,[x(-r),P S (t)}]K(r) 



dT[x,{x{-T)rp s {t)}\K (r) . 



(52) 



Equation ( |52| ) is the ensemble averaged master equation we were looking for. Let us consider the limit of large 
bandwidth and high temperature, 1 < kq C T. Using Eq. (27) we find 



_ r 

K(t) — kq exp 



2tt 



_> r*(r + i|)— >rj(r)+ir|*'(r) 



so that K (r) = r5(r) and if (r) = T(36'{t)/2. We have accordingly, 



dr x(— t) if (r) = — x = Dx 



where we have defined the diffusion coefficient D = Y/2, and 



dr x(—t) K (t) = — 1 — ■ ; 

v ' ' A dr 



= iV^[H s ,x] 

T=Q 4 



IP 



(53) 



(54) 



(55) 



Here we have used [Hs,x] = [p 2 ,x]/2M — —ip/M and defined the friction coefficient 7 — T/3/4M. Collecting 
everything, we finally obtain 



dt 



= [H s ,P S (t)] - D[x, [x,p s (t)}} - i>y[x,{p,p s (t)}] . 



(56) 



The master equation ([56|) is often referred to as the Caldeira-Leggett master equation [[23 24 1. It consists of three 
parts: a von Neumann part which describes the free motion without the coupling to the environment, a diffusive part 
and a dissipative part. Note that the diffusion and friction coefficients satisfy the Einstein relation, D = 2MT^, which 
expresses the fact that diffusion of the particle and damping of its energy have a common physical origin, namely 
the coupling with the heat bath. The Einstein relation is an example of the more general fluctuation-dissipation 
theorem. We also mention that Eq. (|5^) is often derived from an oscillator bath model where the system is coupled 
to a set of independent harmonic oscillators (thus an integrable system). It is quite remarkable that the coupling to a 
complex random matrix environment leads to the same equation. This shows that in the limit considered here (weak 
coupling, high temperature), the master equation (^6|) is independent of the specific structure of the bath as well as 
of the specific form of the coupling and therefore universal |19f| . 



• Coordinate representation and semiclassical limit 

The coordinate representation of the master equation (pq) reads 



dp s {x,x',t) 
dt 



1 



(h s (x) - H s {x'fj - -J (x- xf - ^j(x-x')(p x -p x >) p s (x,x',t) , 



(57) 



where for convenience we have reintroduced the constant h. This corresponds to the master equation derived from 
the influence functional in the limit 



G(x) 



1 ? 
1 - ^x 1 . 



(58) 



The corresponding classical transport equation is obtained by taking the Wigner transform of the density matrix 
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r ipri / r r 
drej q^__j p (g+_, q --,t 



(59) 



Keeping only terms in leading orders of H, this leads to 



dt 



(60) 



This is the usual Klein-Kramers equation. 



C. Relation with the oscillator bath model 

We now turn to a direct comparison of the oscillator bath model and the random matrix model . We still place 
ourselves in the weak coupling limit. The oscillator bath model consists of a particle coupled to a large number of 
independent harmonic oscillators (mass mi and frequency u>i). The Hamiltonian of the composite system is then (see 
p4[ and references therein) 



H = H S 



rriiUJ- 




1 2 N 



m l ujfx 



(61) 



In this model, the coupling Hamiltonian W = x CiXi is bilinear in the position of the system and the the positions 
of the harmonic oscillators. Here Ci are coupling constants. The coupling to the harmonic bath is fully characterized 
by the spectral density function, 



I{w) 



(62) 



In the following we shall relate I(lu) and the variance V a b ■ To this end it is useful to rewrite Eq. (E9J) in the form 



Vab = [p(Sa)p(£b)] 2 £(0J a b = ?b ~ £a) , 



(63) 



where we have introduced the band form factor f(w). It can be shown that with the form ( |63| ) the bath correlation 
function satisfies the KMS condition K(—t) = K(t — i/3) which defines the thermal equilibrium state of the bath. The 
bath correlation function can further be expressed in terms of the spectral density function as 



K(t) 



—I(lu) ^coth ) cos(ujt) ~ isin(o;t) 



(64) 



We see from Eq. (|64|) that the imaginary part K"(t) of the correlation function is simply the Fourier sine transform 
of the spectral density I{uj). Comparing the imaginary parts of Eqs. and (|3), we find, 



/U'l = 2,Tsinh ( ^ ) fU') 



(65) 



This equation provides the desired link between the two models. It shows that there is a one-to-one correspondence 
between the density function I(u>) of the oscillator bath model and the variance V a b 2 of the random band-matrix 
model. It is important to note that both the oscillator bath model and the random matrix model, in the limit of weak 
coupling we consider here, are Gaussian. This means that in these two models the dynamics of the bath is is entirely 
characterized by the two-point correlation function K (t) and it is therefore not necessary to consider higher-order 
correlation functions. 

A straightforward application of Eq. (B5h using the variance Eq. leads to 



= T sinh 







) exp 





(66) 



In the limit lj — > 0, Eq. ( p6[ ) reduces to I{lo) ~ M~/u> which defines the so-called Ohmic regime. In this regime, as we 
already saw in the derivation of the master equation above, the bath correlation function is delta-correlated in time, 
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K (t) = TS(t), which means there are no memory effects. The process in thus Markovian. We can further easily check 
that the band form factor 



2tt 



M a_1 exp 



,2 1 



(67) 



yields I{ui) ~ 7r/3a;f(a;) ~ to a as lu — > 0. This corresponds to fractal (non-Ohmic) environment. Here, contrary to the 
Ohmic case, the bath correlation is given by a inverse power law, 



r (T7T 3 dK' 

K'(t) = -T(a)cos{ — )t- a , K"(t) = \ — 



(68) 



in the limit 1 <C kq <C T. As a result, there are long-time (memory) effects (the process in now non-Markovian) and 
the dynamics of the particle becomes anomalous. We shall develop that point in more detail in the next section. But 
before we would like to show that the band form factor f(w) can be directly related to the velocity autocorrelation 
function (VACF) (v(Q)v(t)) of the system. This is done by using the first fluctuation-dissipation theorem 



C v [z] = (z + -y[z])' 



(69) 



which relates the Laplace transform C v [z] of the normalized VACF C v (t) = (v(0)v(t)) /(v(0) 2 ) to the Laplace transform 
of the damping kernel ^(t) = K (t)/MT (see below). From Eq. (|5l|), we have 



j[z] 



M 



duo e-P"f(aj) 



In the limit z — > 0, the last factor in Eq. (|7^) reduces to a delta function and we obtain 

lim tIzI = — — lim f(z) , 
The final value theorem of the theory of Laplace Transform tells us that 



(70) 



(71) 



lim C v (t) — lim zC v [z] 

t— »oo z—>0 



(72) 



and we can therefore write 



lim C v (t) = lim ( 1 + ^ z-H{z) 



(73) 



We thus see that the long-time behavior of the VACF is determined by the shape of the band form factor of the 
origin. 



D. Langevin equation 

Consider a system operator P (for instance the momentum operator p of the particle). Our aim in the present 
section is to derive a quantum Langevin equation for this operator in the limit of weak coupling. The time evolution 
of the operator P is given by P(t) = exp(iHt) P ex-p(-iHt) (Heisenberg representation). Accordingly, 

P(t) = i[H s (t), P(t)} + i[x(t),P(t)]V(t) . (74) 

This equation is exact. We shall now look for an approximate expression for [x(t), P(t)] = exp(iHt) [x, P] exp(-iHt) 
and V(t) — cxp(iHt)V exjp(-iHt). In the limit of weak coupling, the time evolution operator exp(— iHt) can be 
written in a Dyson series as, 

exp(— iHt) ~ exp(— iH a t) — i / dr cxp(— iH (t — r)) W exp(— iH Q r) . (75) 

Jo 

In lowest order in V, this leads to 
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[x(t),P(t)] = [x(t), P(t)] + i [ dr [x(t - t), [x(t),P(t)}]V(t) , (76) 



" 



and 

r-l 



V(t) = V{t)+ij dTx(t)(v(t-r)V(t) -V(t)V(t-rfj . (77) 
Inserting Eqs. ( |76| ) and ( |77| ) into Eq. ([74]) we then obtain up to second order in V, 

P(t) = i[H s (t),P(t)]- [ dT[x(t-r),[x(t),P(t)]]Rs[V(t-r)V(t)] 



o 



i / dr{x(t-T),[x(t),P(t)]}lm[V(t-T)V(t)] + i[x{t), P(t)]V(t) . (78) 



This equation is valid for any system operator P. We now consider the case where P = p the momentum operator of 
the particle. Using the commutation relation, [x,p] — i, we can rewrite Eq. (Fni) in the form 



p{t) = -U'(x(t)) + 2 [ dr Im[V(t - r)V(t)]x(t - r) - V(t) . (79) 
Jo 

Equation ( f79| ) is almost a Langevin equation. We further note that V(t) depends on the (thermal) initial conditions 
of the heat bath. This makes the force operator £(i) = —V(t) a fluctuating quantity. Introducing first the thermal 
average over the bath and then the average over the random matrix ensemble, we get 



m) = o, mm) - (vwm - m . m 

Moreover, after partial integration, we have 

2 ( dTK"{-T)x(t-T) ~-p [ drK' \t - t)±{t) , (81) 
Jo Jo 

where we have used Eq. (|68|). In the limit of weak coupling and high temperature, the Langevin equation can therefore 
be written as 

MW(t) + M [ drj(t - t)±{t) + U'(x(t)) = £(t) , (82) 
J o 

where we have introduced the damping kernel ^(t) which obeys MT"f(t) = K'(t). This last relation is often referred to 
as the second fluctuation-dissipation theorem. We note that the Langevin equation is completely determined by the 
real part K'{t) of the bath correlation function. Introducing the Riemann-Liouville fractional derivative (— 1<A<0) 

d X f(t) 1 f f(r)dr 

dt> T(-X)J (t-r)^ ' {66 > 

we can rewrite (B2J) in the form of a fractional Langevin equation |E5| 



Mx{t) + M la ^(t) + U'(x(t)) = m , (84) 

where we have defined 7 Q = T(3/(2M sin(a7r/2)) and dropped the overline. The fractional Langevin equation ( |8^ ) 
can be easily solved for the case of a free particle U(x) = by making use of the Laplace transform. We find 

x(t) = x + v B v (t)+ [ dr B v {t - t)£(t) , (85) 
Jo 

where (xq,Vq) are the initial coordinates of the particle and B v (t) = J ' C v (t')dt' is the integral of the (normalized) 
velocity autocorrelation function C v (t). Since j[z] — r y a z a ~ 1 , the Laplace of the VACF is given by 
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By taking the inverse Laplace transform, the velocity autocorrelation function can be written as 



C„(t) = E 2 - a (-l a t 2 - a ) . (87) 
Here we have introduced the Mittag-Leffler function E a (t), which is defined by the series expansion p7| 

00 j-n 

E a {t) = fr, — — TT ■ (88) 

The function E a (t) reduces to the exponential when a = 1. The asymptotic behavior of the Mittag-Leffler function 
( |S8| ) for short and long times is respectively given by ~ exp(f) and ~ — (tr(l — a))^ 1 , < a < 1 and 1 < a < 2. For the 
velocity autocorrelation function (|87]) this yields a typical stretched exponential behavior at short times 

Cv(t) ~ exp ~jf \ , t « , (89) 

r(3-a) ha) L/a 

and an inverse power-law tail at long times 

t a - 2 1 

°v(t) ^7 -r, t > -. —. . (90) 

After time integration, we finally get from Eq. (|87j ) 

B v (t) = tE 2 ^ a . 2 {~ la t 2 ~ a ) , (91) 
where we have used the generalized Mittag-Leffler function E a ^(t) defined as |27]] 

°° f n 
' 1 (an + a) 

In the long-time limit, the generalized Mittag-Leffler function satisfies E a p(t) ~ — (i r(/3— a)) -1 . Accordingly, i? 1 ,(t) 
exhibits a decay of the form 

^Q— 1 

^«(*) ~ — W/ \ 7 when t — > 00 . (93) 
7 Q F(a) 

We emphasize that the solution ( [sB]) of the fractional Langevin equation in the force free case is completely specified 
by the knowledge of the function B v (t) . 

The mean displacement and the mean-square displacement are readily deduced from Eq. (|85|). We find 



(x) =x Q + v tE 2 - a , 2 (- la t 2 ~ a ) ~ (94) 

t^oo Ja T(a) 



a-l 



and 



(x 2 ) = ±Lt 2 E 2 _ a .,{- la t 2 -«) „ Afrn \ z • (95) 
M t^oa 7qM T(l + a) 

In the last equation, thermal initial conditions have been assumed {xq = 0, v 2 = T/M). Equation ( |95| ) shows that 
the coupling to a fractal heat bath leads in general to anomalous diffusion, (x 2 ) ~ t a , a / 1. The band form factor 
( |67| ) gives rise to subdiffusion when a < 1 and to superdiffusion when 1 < a < 2. In general a wide range of anomalous 
transport behaviors can be realized through RMT (quantum chaotic) environments, with the character of the diffusion 
related to the microscopic properties of the quantum environment [20]. 
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VI. RANDOM MATRIX THEORY MASTER EQUATION: 
WEAK COUPLING PLUS TWO BATHS 



We are now in a position to formulate the master equation for the quantum analog of Fig. 1, a system coupled to 
two heat baths B\ and B 2 . The corresponding Hamiltonian is given by 

H = H s + H Bl + H B2 + Q x {x) ®V X + Q 2 {x) <g> V 2 , (96) 

where H$ — J2iPi /^ m i + U{xi,Xj) is a N-particle system Hamiltonian, Qi{x\, otjv) and Q 2 {x\, arjv) are two 
(arbitrary) system operators and V\ and V 2 are two random matrix bath operators. The variance of these two random 
operators is given by Eq. ( p9| ) with the respective parameters of the two baths given by Tfc, Tk and KQ kl k = 1,2. A 
quantum master equation for the system can now be derived using the method described in section VB. This leads 
to the ensemble averaged equation 

^M=- i[Hs ,- Ps (t)} (97) 

poo poo 

dr[Q 1 ,[Qi(~r) 7 p s (t)}}K , 1 (r)- i / dr [Q lt {Qi(-r),p s (t)}} Jf?(r) 



dr[Q 2) [Q 2 (-T),p s (t)}}K 2 (r)-i / dr [Q 2 , {Q 2 (-r), p s (t)}} K 2 (r) 



In the following we take for the system a ID harmonic crystal of unit masses and frequency u>. We attach the two heat 
baths at both ends, X\ and x^, of the linear chain. We thus have Hg — '^2 li pf + uj 2 /2 VJ^ . GijXiXj, where the matrix 
G is defined as Gy = 25^ — — The system operators are respectively given by Qi(xi, ...,xn) = x\ and 

Q 2 [x\, ...,Xff) — xjy. In the Ohmic regime, the corresponding master equation can be easily written in coordinate 
representation in the form 

^ i 1 1 ij 

- L>2 (x N - x' N ) 2 - 72 (x N - x' N ) - I p s (a;, a:', t) , (98) 

with Dk = 2Tfe7fe. Taking the Wigner transform (|5^) we further obtain 

df ( ,v 9/ , ^ at/ 9/ ^ 9 a s 2 / 

2—1 IJ — 1 

+ 2 l2 ^-{ PN f) + 2 72 T 2 |^- . (99) 



Introducing the notation, a;, = lEj+jv — ft, i = 1> N, Eq. ( |99[ ) can be rewritten in the compact form 

r)f 2N r) 1 2N r) 2 



c?t ' 2 ^— ^ dxidx 

i—l — \ 

with £j = VV dijXj. The two 2iV x 2iV matrices a and d are given by 

*=(Jg~r) and d =(S!!)' (101) 

where and I denote the null and unit N x N matrices, and the two N x N matric es R and e obey = 
(271 da + 2 r y 2 SNi)Sij and eg = 2TiRij. The classical generalized Klein-Kramers equation (IOC) has been studied by 
Rieder, Lebowitz and Lieb p8| (the quantum problem has been treated using a Langevin approach in p9[). Of interest 
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here is the stationary non-equilibrium solution df s {xi)/dt = for a small temperature difference, T\ = T + AT/2, 
T 2 =T- AT/2, AT < T. It is given by 



/ s (x) - (2^)- A, Det[b-2]exp 



27V 

-- V 6- 1 



(102) 



where b is the 2iV x 2iV covariance matrix, fry = (xiXj) — J f s (x)xiXjdx.. It is useful to write the covariance matrix 
in the form 



b = 



z' y 



(103) 



where the N x N matrices x, y and z give, respectively, the correlations between the coordinates, the momenta and 
between the coordinates and momenta, — (q-iqj), y-ij — (PiPj), %ij — (q-iPj)- In the limit of large N and small 
coupling 7 = 7i = 72 , one finds 



T 

x=— G"\y 



(T{1) 



V 



T(2) 



T(N) ) 



t ^ 

_7AT 'q 7 AT 



V 



7AT • . 7AT 

a; 2 * m lj 2 

J 



(104) 



Here T(l) = T + AT/2, T(JV) = T - AT/2 and T(i) = T otherwise. We see that the temperature profile in constant 
in the bulk and presents a discontinuity at the edges. By inverting the covariance matrix b we eventually arrive at 



fs(qi,Pi) ~ exp 



E 



pi 



U 7 AT y> p 

9T 9T 2-^, ,,,27 



2T(i) 2T 2T ^w 2 T(i)V%-i %+i 



dU 



dU 



(105) 



This equation is valid in the limit AT, 7 — > 0. We note that Eq. ( |105|) reduces to the equilibrium Boltzmann 
distribution when AT = 0. The last term in (105) is proportional to the heat flux J = 7AT/2 (note that J is not 
proportional to the temperature gradient). Taking the inverse Wigner transform we obtain the stationary density 
operator in the form 



Psilt + = II exp 



1 /7AT 1 
2T(i) V 



2T uj 2 



dU 



'Ji-i 



dU 



H+l 



iT(i)n 



exp 



U 
2T 



(106) 



This form of the solution gives us a guide on how to better solve the general problem of Fig. 1 using a general quantum 
master equation. 



VII. CONCLUSIONS 



The transport properties of classical and quantum systems in non-equilibrium steady states was formulated and 
discussed. While the situation for classical non-equilibrium systems does not seem very well understood yet, it is 
clearly far better understood than the quantum counterpart. At the moment one can formulate a master equation 
for the non-equilibrium density matrix, however, aside from the trivial harmonic case discussed above, very little is 
known. It would be very desirable to see if one can derive Fourier's law or the temperature profile from Random 
Matrix models, and generally how far one can push the matrix models which seem to provide a tractable approach 
to quantum thermalization. 

This work was supported by the Office of Naval Research under contract # N00014-01-1-0594. DK thanks R. 
Olkiewicz and P. Garbaczewki, the organizers of the 38 th Winter School of Theoretical Physics. 
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